R Notebook

here_long <-  -122.3095
here_lat <- 47.6560
seattle = get_map(location = c(here_long, here_lat), zoom = 13, maptype = 'roadmap')

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.656,-122.3095&zoom=13&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false

spd.911 <- spd.911 %>% 
             rowwise() %>% 
             mutate(dist=distVincentyEllipsoid(c(Longitude, Latitude), c(here_long, here_lat)))              

Error in eval(lhs, parent, parent) : object 'spd.911' not found

data <- read.csv('2016-2017-Clean.csv', header = TRUE)
data <- filter(data, !str_detect(Event.Clearance.Description, "HARBOR - DEBRIS, NAVIGATIONAL HAZARDS"))
nrow(data)

[1] 765

ggmap(seattle) +
   geom_point(data = data, aes(x = Longitude, y = Latitude), colour = "red", alpha = 0.75)

  #coord_map()

freq_by_desc <- table(droplevels(data$Event.Clearance.Description))
# View(freq_by_desc)
ggplot(as.data.frame(freq_by_desc), 
       aes(x = Var1, y = Freq)) +
       geom_bar(stat = 'identity') +# create bar plot
    coord_flip()

#Traffic related calls, suspicious circumstances, and disturbances are the the most significant threats to pedestrations
        

ggmap(seattle) +
  geom_point(data = data, aes(x = Longitude, y = Latitude, group = Event.Clearance.Description, color = Event.Clearance.Description), alpha = 0.5, size = 10) +
  facet_wrap(~ Event.Clearance.Description) +
  theme(axis.ticks = element_blank(), 
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        strip.text = element_text(size=50),
        legend.position = "none"
        )

# selecting just ID and location data
df_loc <- data %>% dplyr::select(CAD.CDW.ID, Longitude, Latitude)
# figuring out number of clusters
wss <- c()
# clusters 1 to 15
for (i in 1:15) {
  wss[i] <- sum(kmeans(df_loc, centers=i)$withinss)
}
plot(1:15, wss, type="b", xlab="Number of Clusters",
  ylab="Within groups sum of squares")

# fitting model
fit <- kmeans(df_loc, 10)
fit$centers # look at cluster sizes and means. want clusters to be about equal size

   CAD.CDW.ID Longitude Latitude
1     1966835 -122.3188 47.65906
2     2093969 -122.3142 47.65982
3     2116300 -122.3129 47.66095
4     1899308 -122.3154 47.66025
5     2046956 -122.3160 47.65922
6     2071185 -122.3156 47.66114
7     1864428 -122.3118 47.66149
8     2022152 -122.3172 47.66033
9     1994531 -122.3159 47.66109
10    1932836 -122.3171 47.66133

fit$cluster

  [1]  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
 [40]  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
 [79]  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
[118]  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
[157]  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4 10 10 10 10 10 10 10 10
[196] 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
[235] 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
[274] 10 10 10 10 10 10 10 10 10  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
[313]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
[352]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9
[391]  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9
[430]  9  9  9  9  9  9  9  9  9  9  9  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
[469]  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
[508]  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5
[547]  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6
[586]  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6
[625]  6  6  6  6  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
[664]  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3
[703]  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
[742]  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3

cluster.size <- data.frame(1:10, fit$size)
cluster.size

ggplot(data = cluster.size, aes(x = X1.10, y = fit.size)) +
  geom_bar(stat = 'identity')

ggplot()

ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
aggregate(df_loc, by=list(fit$cluster), FUN=mean)

df_loc

# adding data back into dataframe 
# df_loc <- df_loc %>% mutate(cluster = fit$cluster) 
# View(data)

# hundred block vs TOD
  
by_hr <- table(data$event_clearance_hr)
by_hr

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23 
 38  28  31  52  20  16  17  15  14  16  12  36  20  32  36  39  43  24  25 115  21  38  40  37 

ggplot(as.data.frame(by_hr), aes(x = Var1, y = Freq)) + 
  geom_point() +
  xlab('hour of day')

ggplot(data, aes(x = event_clearance_ts, y = time_until_event_clear)) + 
  geom_point(alpha = 0.25)

ggplot(data, aes(x = Hundred.Block.Location, y = time_until_event_clear)) + 
  geom_point(alpha = 0.25)

  
# selecting just ID and location data
df_loc <- data.w.at.scene %>% dplyr::select(CAD.CDW.ID, Latitude, Longitude)

Error in eval(lhs, parent, parent) : object 'data.w.at.scene' not found

#some useful functions for performing clustering
#extract the lat and long from a dataframe, and run kmeans on it
# x = one of our dataframes
# clusters = how many centers you want kmeans to work with when clustering
fit.clusters <- function(x, clusters) {
  # selecting just ID and location data
  df_loc <- x %>% dplyr::select(CAD.CDW.ID, Latitude, Longitude)
  # fitting model
  fit <- kmeans(df_loc, clusters)
  fit$centers # look at cluster sizes and means. want clusters to be about equal size
  return(fit)
}
#make a plot that will tell you how many clusters might work for a given dataframe
# x = a dataframe
# max = the maximum number of clusters you want to try
find.num.clusters <- function(x, max) {
  if(max > nrow(x)) { stop('Cannot fit more clusters than there are rows in dataframe')}
  df_loc <- x %>% dplyr::select(CAD.CDW.ID, Latitude, Longitude)
  wss = c()
  for (i in 1:max) {
    wss[i] <- sum(kmeans(df_loc, centers=i)$withinss)
  }
  plot(1:max, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")
}
#plot the number of observations in each cluster
# x = a fit object returned from kmeans() or the fit.clusters() function above
plot.cluster.sizes <- function(x) {
  cluster.size <- data.frame(data.frame('clusters' = 1:nrow(x$centers), x$size))
  ggplot(data = cluster.size, aes(x = clusters, y = x.size)) +
  geom_bar(stat = 'identity')
}

Clustering by time of day

morning <- filter(data, 6 <= at_scene_time_hr, at_scene_time_hr < 10 )
mid.day <-  filter(data, 10 <= at_scene_time_hr, at_scene_time_hr < 14 )
afternoon <-  filter(data, 14 <= at_scene_time_hr, at_scene_time_hr < 18 )
evening <-  filter(data, 18 <= at_scene_time_hr, at_scene_time_hr < 22 )
night <-  filter(data, 22 <= at_scene_time_hr | at_scene_time_hr < 2 )
early.morning <-  filter(data, 2 <= at_scene_time_hr, at_scene_time_hr < 6 )
ggmap(seattle) +
  geom_point(data = morning, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = mid.day, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = afternoon, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = evening, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = night, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = early.morning, aes(x = Longitude, y = Latitude), alpha = 0.5)

lengths <- c(nrow(morning), nrow(mid.day), nrow(afternoon), nrow(evening), nrow(night), nrow(early.morning))
names <- c('Morning\n6:00 - 9:59', 'Mid-day\n10:00 - 1:59', 'Afternoon\n2:00 - 5:59', 'Evening\n6:00 - 9:59', 'Night\n10:00 - 1:59', 'Early Morning\n2:00 - 5:59')
by.tod <- data.frame('TOD' = names, 'Count.Crimes' = lengths)
by.tod$TOD = factor(by.tod$TOD, levels = by.tod$TOD)
ggplot(by.tod, aes(x = TOD, y = Count.Crimes)) +
  geom_histogram(stat = 'identity')

Ignoring unknown parameters: binwidth, bins, pad

# find the mode of numeric/character data
Mode <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux)); ux[tab == max(tab)]
}
tod.mean <- mean(data$at_scene_time_hr)
tod.med <- median(data$at_scene_time_hr)
tod.mean

[1] 13.28889

tod.med

[1] 14

Mode(data$at_scene_time_hr)

[1] 17

#What is the most common crime committed at each period?
Mode(morning$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(mid.day$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(afternoon$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(evening$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(night$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(early.morning$Event.Clearance.Description)

[1] CRISIS COMPLAINT - GENERAL
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

#fit kmeans clustering to each time period.
nrow(morning)

[1] 74

find.num.clusters(morning, 10)

fit <- fit.clusters(morning, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(afternoon, 10)

fit <- fit.clusters(mid.day, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

nrow(afternoon)

[1] 176

find.num.clusters(afternoon, 10)

fit <- fit.clusters(afternoon, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(evening, 10)

fit <- fit.clusters(evening, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(night, 10)

fit <- fit.clusters(night, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

#take out general crisis complaint - general
data <- filter(data, Event.Clearance.Description != 'CRISIS COMPLAINT - GENERAL')
morning <- filter(data, 6 <= at_scene_time_hr, at_scene_time_hr < 10 )
mid.day <-  filter(data, 10 <= at_scene_time_hr, at_scene_time_hr < 14 )
afternoon <-  filter(data, 14 <= at_scene_time_hr, at_scene_time_hr < 18 )
evening <-  filter(data, 18 <= at_scene_time_hr, at_scene_time_hr < 22 )
night <-  filter(data, 22 <= at_scene_time_hr | at_scene_time_hr < 2 )
early.morning <-  filter(data, 2 <= at_scene_time_hr, at_scene_time_hr < 6 )
ggmap(seattle) +
  geom_point(data = morning, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = mid.day, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = afternoon, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = evening, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = night, aes(x = Longitude, y = Latitude), alpha = 0.5)

ggmap(seattle) +
  geom_point(data = early.morning, aes(x = Longitude, y = Latitude), alpha = 0.5)

lengths <- c(nrow(morning), nrow(mid.day), nrow(afternoon), nrow(evening), nrow(night), nrow(early.morning))
names <- c('Morning\n6:00 - 9:59', 'Mid-day\n10:00 - 1:59', 'Afternoon\n2:00 - 5:59', 'Evening\n6:00 - 9:59', 'Night\n10:00 - 1:59', 'Early Morning\n2:00 - 5:59')
by.tod <- data.frame('TOD' = names, 'Count.Crimes' = lengths)
by.tod$TOD = factor(by.tod$TOD, levels = by.tod$TOD)
ggplot(by.tod, aes(x = TOD, y = Count.Crimes)) +
  geom_histogram(stat = 'identity')

Ignoring unknown parameters: binwidth, bins, pad

# find the mode of numeric/character data
Mode <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux)); ux[tab == max(tab)]
}
tod.mean <- mean(data$at_scene_time_hr)
tod.med <- median(data$at_scene_time_hr)
tod.mean

[1] 13.69318

tod.med

[1] 14

Mode(data$at_scene_time_hr)

[1] 17

#What is the most common crime committed at each period?
Mode(morning$Event.Clearance.Description)

[1] HAZARDS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(mid.day$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(afternoon$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(evening$Event.Clearance.Description)

[1] HARASSMENT, THREATS
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(night$Event.Clearance.Description)

[1] HARASSMENT, THREATS FIGHT DISTURBANCE  
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

Mode(early.morning$Event.Clearance.Description)

[1] STRONG ARM ROBBERY
7 Levels: ARMED ROBBERY CRISIS COMPLAINT - GENERAL FIGHT DISTURBANCE ... STRONG ARM ROBBERY

#fit kmeans clustering to each time period.
nrow(morning)

[1] 48

find.num.clusters(morning, 10)

fit <- fit.clusters(morning, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(afternoon, 10)

fit <- fit.clusters(mid.day, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

nrow(afternoon)

[1] 107

find.num.clusters(afternoon, 10)

fit <- fit.clusters(afternoon, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(evening, 10)

fit <- fit.clusters(evening, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)

find.num.clusters(night, 10)

fit <- fit.clusters(night, 10)
ggmap(seattle) +
  geom_point(data = as.data.frame(fit$centers), aes(x = Longitude, y = Latitude), alpha = 0.5)

# looking at cluster means
plot.cluster.sizes(fit)